import math

# 地球半径，单位：千米
R = 6371

# 经纬度转换为弧度
def deg_to_rad(deg):
    return deg * (math.pi / 180)

# 计算球面距离
def spherical_distance(lat1, lon1, lat2, lon2):
    lat1_rad = deg_to_rad(lat1)
    lon1_rad = deg_to_rad(lon1)
    lat2_rad = deg_to_rad(lat2)
    lon2_rad = deg_to_rad(lon2)
    d = math.acos(math.sin(lat1_rad) * math.sin(lat2_rad) +
                  math.cos(lat1_rad) * math.cos(lat2_rad) *
                  math.cos(lon1_rad - lon2_rad))
    return R * d

# 经纬度坐标（十进制形式）
points = [(38.7467, 77.6864), (38.745, 77.6933), (38.7422, 77.6911), (38.7289, 77.7103)]

# 计算边长
d12 = spherical_distance(points[0][0], points[0][1], points[1][0], points[1][1])
d23 = spherical_distance(points[1][0], points[1][1], points[2][0], points[2][1])
d34 = spherical_distance(points[2][0], points[2][1], points[3][0], points[3][1])
d41 = spherical_distance(points[3][0], points[3][1], points[0][0], points[0][1])
d13 = spherical_distance(points[0][0], points[0][1], points[2][0], points[2][1])

# 计算三角形面积（海伦公式）
def heron_area(a, b, c):
    p = (a + b + c) / 2
    return math.sqrt(p * (p - a) * (p - b) * (p - c))

# 分割成两个三角形计算面积
S123 = heron_area(d12, d23, d13)
S134 = heron_area(d13, d34, d41)
total_area = S123 + S134
print(f"四个点连线包围的面积约为 {total_area} 平方千米")